
        LOGICAL FUNCTION ptkf_equal(one, two)
C /*
C ** \parambegin
C ** \param{REAL}{one}{floating point number}{IN}
C ** \param{REAL}{two}{floating point number}{IN}
C ** \paramend
C ** \blurb{This function returns TRUE if \pardesc{one} and \pardesc{two}
C ** are equal, or their difference is less than the global 
C ** constant tolerance \pardesc{ptkcpceps}.}
C */
        REAL one, two
        REAL*8 dpone, dptwo
        BYTE ans
        LOGICAL *1 ptk_equal
        external ptk_equal !$PRAGMA C(ptk_equal)

        dpone = one
        dptwo = two
        ans = ptk_equal(%val(dpone), %val(dptwo))
        if (ans .eq. 1) then
                ptkf_equal = .TRUE.
        else
                ptkf_equal = .FALSE.
        endif

        RETURN
        END

        SUBROUTINE ptkf_point(x, y, pt)
C /*
C ** \parambegin
C ** \param{REAL}{x}{x coordinate}{IN}
C ** \param{REAL}{y}{y coordinate}{IN}
C ** \param{REAL}{pt(2)}{real array}{OUT}
C ** \paramend
C ** \blurb{This function puts the values \pardesc{(x,y)} in the
C ** array {\tt pt}.}
C */
        REAL x, y, pt(2)
        pt(1) = x
        pt(2) = y

        RETURN
        END
        
        SUBROUTINE ptkf_point3(x, y, z, pt)
C /*
C ** \parambegin
C ** \param{REAL}{x}{x coordinate}{IN}
C ** \param{REAL}{y}{y coordinate}{IN}
C ** \param{REAL}{z}{z coordinate}{IN}
C ** \param{REAL}{pt(3)}{real array}{OUT}
C ** \paramend
C ** \blurb{This function puts the values \pardesc{(x,y,z)} in the
C ** array {\tt pt}.}
C */
        REAL x, y, z, pt(3)
        
        pt(1) = x
        pt(2) = y
        pt(3) = z

        RETURN
        END
        
        SUBROUTINE ptkf_limit(xmin, xmax, ymin, ymax, lt)
C /*
C ** \parambegin
C ** \param{REAL}{xmin}{minimum x coordinate}{IN}
C ** \param{REAL}{xmax}{maximum x coordinate}{IN}
C ** \param{REAL}{ymin}{minimum y coordinate}{IN}
C ** \param{REAL}{ymax}{maximum y coordinate}{IN}
C ** \param{REAL}{lt(4)}{real array}{OUT}
C ** \paramend
C ** \blurb{This function puts the values \pardesc{(xmin,xmax,ymin,ymax)} 
C ** in the array {\tt lt}.}
C */
        REAL xmin, xmax, ymin, ymax, lt(4)
        
        lt(1) = xmin
        lt(2) = xmax
        lt(3) = ymin
        lt(4) = ymax

        RETURN
        END
        
        SUBROUTINE ptkf_limit3(xmin, xmax, ymin, ymax, zmin, zmax, lt)
C /*
C ** \parambegin
C ** \param{REAL}{xmin}{minimum x coordinate}{IN}
C ** \param{REAL}{xmax}{maximum x coordinate}{IN}
C ** \param{REAL}{ymin}{minimum y coordinate}{IN}
C ** \param{REAL}{ymax}{maximum y coordinate}{IN}
C ** \param{REAL}{zmin}{minimum z coordinate}{IN}
C ** \param{REAL}{zmax}{maximum z coordinate}{IN}
C ** \param{REAL}{lt(6)}{real array}{OUT}
C ** \paramend
C ** \blurb{This function puts the values 
C ** \pardesc{(xmin,xmax,ymin,ymax,zmin,zmax)} 
C ** in the array {\tt lt}.}
C */
        REAL xmin, xmax, ymin, ymax, zmin, zmax, lt(6)
        
        lt(1) = xmin
        lt(2) = xmax
        lt(3) = ymin
        lt(4) = ymax
        lt(5) = zmin
        lt(6) = zmax

        RETURN
        END

        REAL FUNCTION ptkf_dotv3(v1, v2)
C /*
C ** \parambegin
C ** \param{REAL}{v1(3)}{3D vector}{IN}
C ** \param{REAL}{v2(3)}{3D vector}{IN}
C ** \paramend
C ** \blurb{This function evaluates the dot product of the
C **  two 3D vectors \pardesc{v1} and 
C ** \pardesc{v2}, returning it as the value of the function.}
C */
        REAL v1(3), v2(3)
        REAL ptk_dotv3
        external ptk_dotv3 !$PRAGMA C(ptk_dotv3)

        ptkf_dotv3 = ptk_dotv3(v1, v2)

        RETURN
        END

        REAL FUNCTION ptkf_dotv(v1, v2)
C /*
C ** \parambegin
C ** \param{REAL}{v1(2)}{2D vector}{IN}
C ** \param{REAL}{v2(2)}{2D vector}{IN}
C ** \paramend
C ** \blurb{Evaluates the dot product of the two 2D vectors \pardesc{v1} and
C ** \pardesc{v2}, returning it as the value of the function.}
C */
        REAL v1(2), v2(2)
        REAL ptk_dotv
        external ptk_dotv !$PRAGMA C(ptk_dotv)
      
        ptkf_dotv = ptk_dotv(v1, v2)

        RETURN
        END

        SUBROUTINE ptkf_crossv3(v1, v2, v3)
C /*
C ** \parambegin
C ** \param{REAL}{v1(3)}{3D vector}{IN}
C ** \param{REAL}{v2(3)}{3D vector}{IN}
C ** \param{REAL}{v3(3)}{3D vector}{OUT}
C ** \paramend
C ** \blurb{This function evaluates the cross product of
C ** the two 3D vectors \pardesc{v1} and
C ** \pardesc{v2}, returning the new vector as
C ** the function result. Since a local copy is made statements such as
C ** {\tt call ptkf\_crossv(v1, v2, v2)}
C ** will produce the correct answer.}
C */
        REAL v1(3), v2(3), v3(3)
        REAL temp(3)         

        temp(1) = v1(2) * v2(3) - v1(3) * v2(2)
        temp(2) = v1(3) * v2(1) - v1(1) * v2(3)
        temp(3) = v1(1) * v2(2) - v1(2) * v2(1)
        v3(1) = temp(1)
        v3(2) = temp(2)
        v3(3) = temp(3)

        RETURN
        END
                
        LOGICAL FUNCTION ptkf_nullv3(vec)
C /*
C ** \parambegin
C ** \param{REAL}{vec(3)}{3D vector}{IN}
C ** \paramend
C ** \blurb{This function returns \pardesc{TRUE} if the modulus
C ** of the 3D vector\pardesc{vec} 
C ** is less than the global tolerance \pardesc{ptkpceps}, otherwise 
C ** \pardesc{FALSE}.}
C */
        REAL vec(3)
        BYTE ans
        LOGICAL *1 ptk_nullv3
        external ptk_nullv3 !$PRAGMA C(ptk_nullv3)

        ans = ptk_nullv3(vec)
        if (ans .eq. 1) then
                ptkf_nullv3 = .TRUE.
        else
                ptkf_nullv3 = .FALSE.
        endif
        
        RETURN
        END
        
        LOGICAL FUNCTION ptkf_nullv(vec)
C /*
C ** \parambegin
C ** \param{REAL}{vec(2)}{2D vector}{IN}
C ** \paramend
C ** \blurb{This function returns \pardesc{TRUE} if the
C ** modulus of the 2D vector \pardesc{vec} 
C ** is less than the global tolerance \pardesc{ptkpceps}, otherwise 
C ** \pardesc{FALSE}.}
C */
        REAL vec(3)
        BYTE ans
        LOGICAL *1 ptk_nullv
        external ptk_nullv !$PRAGMA C(ptk_nullv)

        ans = ptk_nullv(vec)
        if (ans .eq. 1) then
                ptkf_nullv = .TRUE.
        else
                ptkf_nullv = .FALSE.
        endif
        
        RETURN
        END
        
        REAL FUNCTION ptkf_modv3(vec)
C /*
C ** \parambegin
C ** \param{REAL}{vec(3)}{3D vector}{IN}
C ** \paramend
C ** \blurb{Returns the modulus of the vector \pardesc{vec}.}
C */
        REAL vec(3)
        REAL ptk_modv3
        external ptk_modv3 !$PRAGMA C(ptk_modv3)

        ptkf_modv3 = ptk_modv3(vec)
        
        RETURN
        END
        
        REAL FUNCTION ptkf_modv(vec)
C /*
C ** \parambegin
C ** \param{REAL}{vec(2)}{2D vector}{IN}
C ** \paramend
C ** \blurb{This function returns the modulus of the 2D vector \pardesc{vec}.}
C */
        REAL vec(2)
        REAL ptk_modv
        external ptk_modv !$PRAGMA C(ptk_modv)
        
        ptkf_modv = ptk_modv(vec)
        
        RETURN
        END
        
        SUBROUTINE ptkf_unitv3(vec, uvec)
C /*
C ** \parambegin
C ** \param{REAL}{vec(3)}{3D vector}{IN}
C ** \param{REAL}{uvec(3)}{3D vector}{OUT}
C ** \paramend
C ** \blurb{This function generates and returns a unit
C ** vector in {\tt uvec} from the supplied 3D vector 
C ** \pardesc{vec}.}
C */
        REAL vec(3), uvec(3), modu
        REAL ptkf_modv3
        LOGICAL ptkf_equal

        modu = ptkf_modv3(vec)
        if (ptkf_equal(modu, 0.0) .eq. .FALSE.) then
          call ptkf_point3(vec(1) / modu, vec(2) / modu, vec(3) / modu, 
& uvec)
        else 
          call ptkf_point3(0.0, 0.0, 0.0, uvec)
        endif
        
        RETURN
        END
        
        SUBROUTINE ptkf_unitv(vec, uvec)
C /*
C ** \parambegin
C ** \param{REAL}{vec(2)}{2D vector}{IN}
C ** \param{REAL}{uvec(2)}{2D vector}{OUT}
C ** \paramend
C ** \blurb{This function generates and returns a unit vector 
C ** in {\tt uvec} from the supplied 2D vector \pardesc{vec}.}
C */
        REAL vec(2), uvec(2), modu
        REAL ptkf_modv
        LOGICAL ptkf_equal

        modu = ptkf_modv(vec)
        if (ptkf_equal(modu, 0.0) .eq. .FALSE.) then
          call ptkf_point(vec(1) / modu, vec(2) / modu, uvec)
        else 
          call ptkf_point(0.0, 0.0, uvec)
        endif        
        
        RETURN
        END
        
        SUBROUTINE ptkf_scalev3(vec, scale, svec)
C /*
C ** \parambegin
C ** \param{REAL}{vec(3)}{3D vector}{IN}
C ** \param{REAL}{scale}{scale factor}{IN}
C ** \param{REAL}{svec(3)}{3D vector}{OUT}
C ** \paramend
C ** \blurb{This function multiplies the 3D vector \pardesc{vec}
C **  by the  scalar \pardesc{scale} and
C ** returns the result in {\tt svec}.}
C */
        REAL vec(3), scale, svec(3)

        call ptkf_point3(vec(1) * scale, vec(2) * scale, vec(3) * scale, 
& svec)
       
        RETURN
        END
        
        SUBROUTINE ptkf_scalev(vec, scale, svec)
C /*
C ** \parambegin
C ** \param{REAL}{vec(2)}{2D vector}{IN}
C ** \param{REAL}{scale}{scale factor}{IN}
C ** \param{REAL}{svec(2)}{2D vector}{OUT}
C ** \paramend
C ** \blurb{This function multiplies the 2D vector
C ** \pardesc{vec} by the scalar \pardesc{scale} and
C ** returns the result in {\tt svec}.}
C */
        REAL vec(2), scale, svec(2)

        call ptkf_point(vec(1) * scale, vec(2) * scale, svec)

        RETURN
        END
        
        SUBROUTINE ptkf_subv3(p1, p2, p3)
C /*
C ** \parambegin
C ** \param{REAL}{p1(3)}{3D vector}{IN}
C ** \param{REAL}{p2(3)}{3D vector}{IN}
C ** \param{REAL}{p3(3)}{3D vector}{OUT}
C ** \paramend
C ** \blurb{This function evaluates the 3D vector \pardesc {p1-p2} and
C ** returns the result in {\tt p3}.}
C */
        REAL p1(3), p2(3), p3(3)
        
        call ptkf_point3(p1(1) - p2(1), p1(2)- p2(2), p1(3) - p2(3), 
& p3)
        
        RETURN
        END
        
        SUBROUTINE ptkf_subv(p1, p2, p3)
C /*
C ** \parambegin
C ** \param{REAL}{p1(2)}{2D vector}{IN}
C ** \param{REAL}{p2(2)}{2D vector}{IN}
C ** \param{REAL}{p3(2)}{2D vector}{OUT}
C ** \paramend
C ** \blurb{This function evaluates the 2D vector \pardesc {p1-p2}
C ** and returns the result in {\tt p3}.}
C */
        REAL p1(2), p2(2), p3(2)
                
        call ptkf_point(p1(1) - p2(1), p1(2)- p2(2), p3)

        RETURN
        END
        
        SUBROUTINE ptkf_addv3(p1, p2, p3)
C /*
C ** \parambegin
C ** \param{REAL}{p1(3)}{3D vector}{IN}
C ** \param{REAL}{p2(3)}{3D vector}{IN}
C ** \param{REAL}{p3(3)}{3D vector}{OUT}
C ** \paramend
C ** \blurb{This function evaluates the 3D vector \pardesc{p1+p3} and 
C ** returns the result in {\tt p3}.}
C */
        REAL p1(3), p2(3), p3(3)
        
        call ptkf_point3(p1(1) + p2(1), p1(2) + p2(2), p1(3) + p2(3), 
& p3)
        
        RETURN
        END
        
        SUBROUTINE ptkf_addv(p1, p2, p3)
C /*
C ** \parambegin
C ** \param{REAL}{p1(2)}{2D vector}{IN}
C ** \param{REAL}{p2(2)}{2D vector}{IN}
C ** \param{REAL}{p3(2)}{2D vector}{OUT}
C ** \paramend
C ** \blurb{This function evaluates the 2D vector \pardesc{p1+p2}
C ** and returns the result in {\tt p3}.}
C */
        REAL p1(2), p2(2), p3(2)
        
        call ptkf_point(p1(1) + p2(1), p1(2) + p2(2), p3)
        
        RETURN
        END
        
        SUBROUTINE ptkf_unitmatrix(matrix)
C /*
C ** \parambegin
C ** \param{REAL}{matrix(3, 3)}{3x3 matrix}{IN}
C ** \paramend
C ** \blurb{This procedure creates a unit $3\times 3$ matrix, and stores 
C ** it in \pardesc{matrix}.}
C */
        REAL matrix(3,3)
        external ptk_unitmatrix !$PRAGMA C(ptk_unitmatrix)
        
        call  ptk_unitmatrix(matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_unitmatrix3(matrix)
C /*
C ** \parambegin
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{IN}
C ** \paramend
C ** \blurb{This procedure creates a unit $4\times 4$ matrix, and stores 
C ** it in \pardesc{matrix}.}
C */
        REAL matrix(4,4)
        external ptk_unitmatrix3 !$PRAGMA C(ptk_unitmatrix3)
        
        call  ptk_unitmatrix3 (matrix)
        
        RETURN
        END
                
        SUBROUTINE ptkf_transposematrix3(matrix, result)
C /*
C ** \parambegin
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{IN}
C ** \param{REAL}{result(4, 4)}{4x4 matrix}{OUT}
C ** \paramend
C ** \blurb{This function transposes \pardesc{matrix}, and returns the result
C ** in \pardesc{result}.
C ** Note that \pardesc{result} can be the same variable
C ** as \pardesc{matrix} since a copy is made 
C ** first.}
C */
        REAL matrix(4,4), result(4,4)
        external ptk_transposematrix3 !$PRAGMA C(ptk_transposematrix3)
        
        call  ptk_transposematrix3(matrix, result)
        
        RETURN
        END
        
        SUBROUTINE ptkf_transposematrix(matrix, result)
C /*
C ** \parambegin
C ** \param{REAL}{matrix(3, 3)}{3x3 matrix}{IN}
C ** \param{REAL}{result(3, 3)}{3x3 matrix}{OUT}
C ** \paramend
C ** \blurb{This function transposes \pardesc{matrix}, and returns the result
C ** in \pardesc{result}.
C ** Note that \pardesc{result} can be the same variable
C ** as \pardesc{matrix} since a copy is made 
C ** first.}
C */
        REAL matrix(3,3), result(3,3)
        external ptk_transposematrix !$PRAGMA C(ptk_transposematrix)

        call  ptk_transposematrix(matrix, result)
        
        RETURN
        END
        
        SUBROUTINE ptkf_multiplymatrix3(matrix1, matrix2, result)
C /*
C ** \parambegin
C ** \param{REAL}{matrix1(4, 4)}{4x4 matrix}{IN}
C ** \param{REAL}{matrix2(4, 4)}{4x4 matrix}{IN}
C ** \param{REAL}{result(4, 4)}{4x4 matrix}{OUT}
C ** \paramend
C ** \blurb{This function makes \pardesc{result} the product of
C ** the $4 \times 4$ matrices  \pardesc{matrix1} and
C ** \pardesc{matrix2}, with {\tt result $\leftarrow$ matrix1 * matrix2}.
C ** Note that \pardesc{result} can also be \pardesc{matrix1} or 
C ** \pardesc{matrix2} since a copy is made.}
C */
        REAL matrix1(4,4), matrix2(4,4), result(4,4)
        external ptk_multiplymatrix3 !$PRAGMA C(ptk_multiplymatrix3)

        call  ptk_multiplymatrix3(matrix1, matrix2, result)
        
        RETURN
        END
        
        SUBROUTINE ptkf_multiplymatrix(matrix1, matrix2, result)
C /*
C ** \parambegin
C ** \param{REAL}{matrix1(3, 3)}{3x3 matrix}{IN}
C ** \param{REAL}{matrix2(3, 3)}{3x3 matrix}{IN}
C ** \param{REAL}{result(3, 3)}{3x3 matrix}{OUT}
C ** \paramend
C ** \blurb{This function makes \pardesc{result} the product of the
C ** $3 \times 3$ matrices \pardesc{matrix1} and
C ** \pardesc{matrix2}, with {\tt result $\leftarrow$ matrix1 * matrix2}.
C ** Note that \pardesc{result} can also be \pardesc{matrix1} or 
C ** \pardesc{matrix2} since a copy is made.}
C */
        REAL matrix1(3,3), matrix2(3,3), result(3,3)
        external ptk_multiplymatrix !$PRAGMA C(ptk_multiplymatrix)

        call  ptk_multiplymatrix(matrix1, matrix2, result)
        
        RETURN
        END
        
        SUBROUTINE ptkf_concatenatematrix3(operation, matrix1, matrix2, 
& result)
C /*
C ** \parambegin
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix1(4, 4)}{4x4 matrix}{IN}
C ** \param{REAL}{matrix2(4, 4)}{4x4 matrix}{IN}
C ** \param{REAL}{result(4, 4)}{4x4 matrix}{OUT}
C ** \paramend
C ** \blurb{This function concatenates the $4 \times 4$ matrices
C ** \pardesc{matrix1}  and \pardesc{matrix2}
C ** on the basis of \pardesc{operation}.
C ** The result is stored in \pardesc{result}.
C ** Note that \pardesc{result} can also be \pardesc{matrix1} or 
C ** \pardesc{matrix2}
C ** since a copy is made. When \pardesc{operation} is 
C ** \pardesc{preconcatenate}, then 
C ** \pardesc{result $\leftarrow$ matrix2 * matrix1}.
C ** When \pardesc{operation} is \pardesc{postconcatenate}, 
C ** \pardesc{result $\leftarrow$  matrix1 * matrix2}.}
C */
        INTEGER operation
        REAL matrix1(4,4), matrix2(4,4), result(4,4)
        external ptk_concatenatematrix3 
& !$PRAGMA C(ptk_concatenatematrix3)

        call  ptk_concatenatematrix3(%val(operation), matrix1, matrix2, 
& result)
         
        RETURN
        END
   
        SUBROUTINE ptkf_concatenatematrix(operation, matrix1, matrix2, 
& result)
C /*
C ** \parambegin
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix1(3, 3)}{3x3 matrix}{IN}
C ** \param{REAL}{matrix2(3, 3)}{3x3 matrix}{IN}
C ** \param{REAL}{result(3, 3)}{3x3 matrix}{OUT}
C ** \paramend
C ** \blurb{This function concatenates
C ** the $3 \times 3$ matrices \pardesc{matrix1}  and \pardesc{matrix2}
C ** on the basis of \pardesc{operation}.
C ** The result is stored in \pardesc{result}.
C ** Note that \pardesc{result} can also be \pardesc{matrix1} or 
C ** \pardesc{matrix2}
C ** since a copy is made. When \pardesc{operation} is
C ** \pardesc{preconcatenate}, then 
C ** \pardesc{result $\leftarrow$ matrix2 * matrix1}.
C ** When \pardesc{operation} is \pardesc{postconcatenate}, 
C ** \pardesc{result $\leftarrow$ matrix1 * matrix2}.}
C */
        INTEGER operation
        REAL matrix1(3,3), matrix2(3,3), result(3,3)
        external ptk_concatenatematrix !$PRAGMA C(ptk_concatenatematrix)
        
        call  ptk_concatenatematrix(%val(operation), matrix1, matrix2, 
& result)
        
        RETURN
        END
        
        SUBROUTINE ptkf_shift3(shift, operation, matrix)
C /*
C ** \parambegin
C ** \param{REAL}{shift(3)}{shift vector}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{OUT}
C ** \paramend
C ** \blurb{This function computes
C ** a matrix to perform the specified 3D shift and concatenates 
C ** this matrix with \pardesc{matrix} on the basis of \pardesc{operation}.}
C */
        REAL shift(3)
        INTEGER operation
        REAL matrix(4,4)
        external ptk_shift3 !$PRAGMA C(ptk_shift3)
        
        call  ptk_shift3(shift, %val(operation), matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_shift(shift, operation, matrix)
C /*
C ** \parambegin
C ** \param{REAL}{shift(2)}{shift vector}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(3, 3)}{3x3 matrix}{OUT}
C ** \paramend
C ** \blurb{This function computes
C ** a matrix to perform the specified 2D shift and concatenates 
C ** this matrix with \pardesc{matrix} on the basis of \pardesc{operation}.}
C */
        REAL shift(2)
        INTEGER operation
        REAL matrix(3,3)
        external ptk_shift !$PRAGMA C(ptk_shift)
        
        call  ptk_shift(shift, %val(operation), matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_scale3(scale, operation, matrix)
C /*
C ** \parambegin
C ** \param{REAL}{scale(3)}{scale vector}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{OUT}
C ** \paramend
C ** \blurb{This function 
C ** computes a matrix to perform the specified 3D scale and concatenates 
C ** this with \pardesc{matrix} on the basis of \pardesc{operation}.}
C */
        REAL scale(3)
        INTEGER operation
        REAL matrix(4,4)
        external ptk_scale3 !$PRAGMA C(ptk_scale3)
        
        call  ptk_scale3(scale, %val(operation), matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_scale(scale, operation, matrix)
C /*
C ** \parambegin
C ** \param{REAL}{scale(2)}{scale vector}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(3, 3)}{3x3 matrix}{OUT}
C ** \paramend
C ** \blurb{This function computes
C **  a matrix to perform the specified 2D scale and concatenates 
C ** this matrix with \pardesc{matrix} on the basis of \pardesc{operation}.}
C */
        REAL scale(2)
        INTEGER operation
        REAL matrix(3,3)
        external ptk_scale !$PRAGMA C(ptk_scale)
        
        call  ptk_scale(scale, %val(operation), matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_rotatecs3(costheta, sinetheta, axis, operation, 
& matrix)
C /*
C ** \parambegin
C ** \param{REAL}{costheta}{cosine of angle}{IN}
C ** \param{REAL}{sinetheta}{sine of angle}{IN}
C ** \param{INTEGER}{axis}{x, y or z axis}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{OUT}
C ** \paramend
C ** \blurb{This function computes
C ** a matrix to perform the specified 3D rotation and concatenates 
C ** this matrix with 
C ** \pardesc{matrix} on the basis of \pardesc{operation}.
C ** This form assumes that the rotation is specified
C ** using the $\cos(theta)$ and $\sin(theta)$ terms.
C ** Note that no check is made to ensure that the sum of the squares of these
C ** terms is 1.}
C */
        REAL costheta, sinetheta
        INTEGER axis, operation
        REAL matrix(4,4)
        REAL*8 dpcostheta, dpsinetheta
        external ptk_rotatecs3 !$PRAGMA C(ptk_rotatecs3)
        
        dpcostheta = costheta
        dpsinetheta = sinetheta
        call  ptk_rotatecs3(%val(dpcostheta), %val(dpsinetheta), 
& %val(axis), %val(operation), matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_rotatecs(costheta, sinetheta, axis, operation, 
& matrix)
C /*
C ** \parambegin
C ** \param{REAL}{costheta}{cosine of angle}{IN}
C ** \param{REAL}{sinetheta}{sine of angle}{IN}
C ** \param{INTEGER}{axis}{x, y or z axis}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(3, 3)}{3x3 matrix}{OUT}
C ** \paramend
C ** \blurb{This function computes
C ** a matrix to perform the specified 2D rotation and concatenates 
C ** this matrix with 
C ** \pardesc{matrix} on the basis of \pardesc{operation}.
C ** This form assumes that the rotation is specified
C ** using the $\cos(theta)$ and $\sin(theta)$ terms.
C ** Note that no check is made to ensure that the sum of the squares of these
C ** terms is 1.}
C */
        REAL costheta, sinetheta
        INTEGER axis, operation
        REAL matrix(3,3)
        REAL*8 dpcostheta, dpsinetheta
        external ptk_rotatecs !$PRAGMA C(ptk_rotatecs)
        
        dpcostheta = costheta
        dpsinetheta = sinetheta
        call  ptk_rotatecs(%val(dpcostheta), %val(dpsinetheta), 
& %val(axis), %val(operation), matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_rotate3(rotation, axis, operation, matrix)
C /*
C ** \parambegin
C ** \param{REAL}{rotation}{angle in degrees}{IN}
C ** \param{INTEGER}{axis}{x, y or z axis}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{OUT}
C ** \paramend
C ** \blurb{This function computes
C ** a matrix to perform the specified 3D rotation and concatenates 
C ** this matrix with 
C ** \pardesc{matrix} on the basis of \pardesc{operation}.
C ** \pardesc{rotation} is expressed in degrees.}
C */
        REAL rotation
        INTEGER axis, operation
        REAL matrix(4,4)
        REAL*8 dprotation

        external ptk_rotate3 !$PRAGMA C(ptk_rotate3)
        
        dprotation = rotation
        call  ptk_rotate3(%val(dprotation), %val(axis), 
& %val(operation), matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_rotate(rotation, axis, operation, matrix)
C /*
C ** \parambegin
C ** \param{REAL}{rotation}{angle in degrees}{IN}
C ** \param{INTEGER}{axis}{x, y or z axis}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(3, 3)}{3x3 matrix}{OUT}
C ** \paramend
C ** \blurb{This function computes
C ** a matrix to perform the specified 2D rotation and concatenates 
C ** this matrix with 
C ** \pardesc{matrix} on the basis of \pardesc{operation}.
C ** \pardesc{rotation} is expressed in degrees.}
C */
        REAL rotation
        INTEGER axis, operation
        REAL matrix(3,3)
        REAL*8 dprotation
        external ptk_rotate !$PRAGMA C(ptk_rotate)
        
        dprotation = rotation
        call  ptk_rotate(%val(dprotation), %val(axis), %val(operation), 
& matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_shear3(shearaxis, sheardir, shearfactor, 
& operation, matrix)
C /*
C ** \parambegin
C ** \param{INTEGER}{shearaxis}{x, y or z axis}{IN}
C ** \param{INTEGER}{sheardir}{x, y or z direction}{IN}
C ** \param{REAL}{shearfactor}{shear factor}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{OUT}
C ** \paramend
C ** \blurb{This function computes a matrix to perform the specified 3D
C **  shear and concatenates 
C ** this matrix with 
C ** \pardesc{matrix} on the basis of \pardesc{operation}.
C ** The shear is specified as an amount \pardesc{f} about axis \pardesc{i}
C ** in direction \pardesc{j}.}
C */
        INTEGER shearaxis, sheardir
        REAL shearfactor
        INTEGER operation
        REAL matrix(4,4)
        REAL*8 dpshearfactor
        external ptk_shear3 !$PRAGMA C(ptk_shear3)
        
        dpshearfactor = shearfactor
        call  ptk_shear3(%val(shearaxis), %val(sheardir), 
& %val(dpshearfactor), %val(operation), matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_shear(shearaxis, sheardir, shearfactor, 
& operation, matrix)
C /*
C ** \parambegin
C ** \param{INTEGER}{shearaxis}{x or y axis}{IN}
C ** \param{INTEGER}{sheardir}{x or y direction}{IN}
C ** \param{REAL}{shearfactor}{shear factor}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(3, 3)}{3x3 matrix}{OUT}
C ** \paramend
C ** \blurb{This function computes a matrix to perform the specified 2D
C **  shear and concatenates 
C ** this matrix with 
C ** \pardesc{matrix} on the basis of \pardesc{operation}.
C ** The shear is specified as an amount \pardesc{shearfactor} 
C ** about axis \pardesc{shearaxis}
C ** in direction \pardesc{sheardir}.}
C */
        INTEGER shearaxis, sheardir
        REAL shearfactor
        INTEGER operation
        REAL matrix(3,3)
        REAL*8 dpshearfactor
        external ptk_shear !$PRAGMA C(ptk_shear)
        
        dpshearfactor = shearfactor
        call  ptk_shear(%val(shearaxis), %val(sheardir), 
& %val(dpshearfactor), %val(operation), matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_rotatevv3(v1, v2, operation, matrix, error)
C /*
C ** \parambegin
C ** \param{REAL}{v1(3)}{3D vector}{IN}
C ** \param{REAL}{v2(3)}{3D vector}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{OUT}
C ** \param{INTEGER}{error}{error code}{OUT}
C ** \paramend
C ** \blurb{This function computes
C **  a matrix to perform the rotation (about the origin) of the 3D vector
C ** \pardesc{v1} to the 3D vector 
C ** \pardesc{v2}, and concatenates this matrix with 
C ** \pardesc{matrix} on the basis of \pardesc{operation} (\cite{rog:mecg}, 
C **  pages 55--59). If the parameters are invalid, \pardesc{error} is set to 
C ** -1. Otherwise, its value is \pardesc{ptkcpcok}.}
C */
        REAL v1(3), v2(3)
        INTEGER operation
        REAL matrix(4,4)
        INTEGER error
        external ptk_rotatevv3 !$PRAGMA C(ptk_rotatevv3)
        
        call  ptk_rotatevv3(v1, v2, %val(operation), matrix, error)
        
        RETURN
        END
        
        SUBROUTINE ptkf_rotatevv(v1, v2, operation, matrix, error)
C /*
C ** \parambegin
C ** \param{REAL}{v1(2)}{2D vector}{IN}
C ** \param{REAL}{v2(2)}{2D vector}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(3, 3)}{3x3 matrix}{OUT}
C ** \param{INTEGER}{error}{error code}{OUT}
C ** \paramend
C ** \blurb{This function computes
C **  a matrix to perform the rotation (about the origin) of the 3D vector
C ** \pardesc{v1} to the 3D vector 
C ** \pardesc{v2}, and concatenates this matrix with 
C ** \pardesc{matrix} on the basis of \pardesc{operation} (\cite{rog:mecg}, 
C **  pages 55--59). If the parameters are invalid, \pardesc{error} is set to 
C ** -1. Otherwise, its value is \pardesc{ptkcpcok}.}
C */
        REAL v1(2), v2(2)
        INTEGER operation
        REAL matrix(3,3)
        INTEGER error
        external ptk_rotatevv !$PRAGMA C(ptk_rotatevv)
        
        call  ptk_rotatevv(v1, v2, %val(operation), matrix, error)
        
        RETURN
        END
        
        SUBROUTINE ptkf_rotateline3(p1, p2, theta, operation, matrix, 
& error)
C /*
C ** \parambegin
C ** \param{REAL}{p1(3)}{3D point}{IN}
C ** \param{REAL}{p2(3)}{3D point}{IN}
C ** \param{REAL}{theta}{angle in degrees}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{OUT}
C ** \param{INTEGER}{error}{error code}{OUT}
C ** \paramend
C ** \blurb{This function computes
C ** a matrix to perform a 3D rotation of
C ** \pardesc{theta} degrees
C ** about the line connecting \pardesc{p1} to \pardesc{p2},
C ** and concatenates this matrix with 
C ** \pardesc{matrix} on the basis of \pardesc{operation}.
C ** If the parameters are invalid, \pardesc{error} is set to 
C ** -1. Otherwise, its value is \pardesc{ptkcpcok}.}
C */
        REAL p1(3), p2(3), theta
        INTEGER operation
        REAL matrix(4,4)
        INTEGER error
        REAL*8 dptheta
        external ptk_rotateline3 !$PRAGMA C(ptk_rotateline3)

        dptheta = theta
        call  ptk_rotateline3(p1, p2, %val(dptheta), %val(operation), 
& matrix, error)
        
        RETURN
        END
        
        SUBROUTINE ptkf_rotateline(p1, p2, theta, operation, matrix, 
& error)
C /*
C ** \parambegin
C ** \param{REAL}{p1(2)}{2D point}{IN}
C ** \param{REAL}{p2(2)}{2D point}{IN}
C ** \param{REAL}{theta}{angle in degrees}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(3, 3)}{3x3 matrix}{OUT}
C ** \param{INTEGER}{error}{error code}{OUT}
C ** \paramend
C ** \blurb{This function computes
C ** a matrix to perform a 2D rotation of
C ** \pardesc{theta} degrees
C ** about the line connecting \pardesc{p1} to \pardesc{p2},
C ** and concatenates this matrix with 
C ** \pardesc{matrix} on the basis of \pardesc{operation}.
C ** If the parameters are invalid, \pardesc{error} is set to 
C ** -1. Otherwise, its value is \pardesc{ptkcpcok}.}
C */
        REAL p1(2), p2(2), theta
        INTEGER operation
        REAL matrix(3,3)
        INTEGER error
        REAL*8 dptheta
        external ptk_rotateline !$PRAGMA C(ptk_rotateline)
        
        dptheta = theta
        call  ptk_rotateline(p1, p2, %val(dptheta), %val(operation), 
& matrix, error)
        
        RETURN
        END
        
        SUBROUTINE ptkf_pt3topt4(pt, pt4)
C /*
C ** \parambegin
C ** \param{REAL}{pt(3)}{3D point}{IN}
C ** \param{REAL}{pt4(4)}{4D point}{OUT}
C ** \paramend
C ** \blurb{This function converts the 3D point \pardesc{pt} to a 4D point,
C ** {\tt pt4}. The $w$ coordinate of the 
C ** 4D point is set to $1.0$.}
C */
        REAL pt(3), pt4(4)

        pt4(1) = pt(1)
        pt4(2) = pt(2)
        pt4(3) = pt(3)
        pt4(4) = 1.0
        
        RETURN
        END
        
        SUBROUTINE ptkf_pt4topt3(pt, pt3)
C /*
C ** \parambegin
C ** \param{REAL}{pt(4)}{4D point}{IN}
C ** \param{REAL}{pt3(3)}{3D point}{OUT}
C ** \paramend
C ** \blurb{This function converts the 4D point \pardesc{pt} to a 3D point
C ** {\tt pt3}, by dividing by $w$.}
C */
        REAL pt(4), pt3(3), w
        LOGICAL ptkf_equal, ans
  
        ans = ptkf_equal(pt(4), 0.0)
        if (ans .eq. .TRUE.) then
          w = 1.0 / 1.0e-7
        else
          w = 1.0 / pt(4)
        endif

        pt3(1) = pt(1) * w
        pt3(2) = pt(2) * w
        pt3(3) = pt(3) * w
        
        RETURN
        END
        
        SUBROUTINE ptkf_transform4(matrix, point, tpoint)
C /*
C ** \parambegin
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{IN}
C ** \param{REAL}{point(4)}{4D point}{IN}
C ** \param{REAL}{tpoint(4)}{4D point}{OUT}
C ** \paramend
C ** \blurb{This function performs the 4D transformation 
C ** (that is, with no homogeneous division) of the
C ** point \pardesc{point} by the $4 \times 4$ matrix
C ** \pardesc{matrix}, and returns the result in {\tt tpoint}.}
C */
        REAL matrix(4,4), point(4), tpoint(4)

        tpoint(1) = matrix(1, 1) * point(1) + matrix(1, 2) * point(2) + 
& matrix(1, 3) * point(3) + matrix(1, 4) * point(4)
        tpoint(2) = matrix(2, 1) * point(1) + matrix(2, 2) * point(2) + 
& matrix(2, 3) * point(3) + matrix(2, 4) * point(4)
        tpoint(3) = matrix(3, 1) * point(1) + matrix(3, 2) * point(2) + 
& matrix(3, 3) * point(3) + matrix(3, 4) * point(4)
        tpoint(4) = matrix(4, 1) * point(1) + matrix(4, 2) * point(2) + 
& matrix(4, 3) * point(3) + matrix(4, 4) * point(4)
        
        RETURN
        END
        
        SUBROUTINE ptkf_transform3(matrix, point, tpoint)
C /*
C ** \parambegin
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{IN}
C ** \param{REAL}{point(3)}{3D point}{IN}
C ** \param{REAL}{tpoint(3)}{3D point}{OUT}
C ** \paramend
C ** \blurb{This function transforms the 3D point \pardesc{point} by
C ** the $4 \times 4$ matrix \pardesc{matrix},
C ** including homogeneous division. 
C ** The result is returned in {\tt tpoint}.}
C */
        REAL matrix(4,4), point(3), tpoint(3), temp(4)

        temp(1) = matrix(1, 1) * point(1) + matrix(1, 2) * point(2) + 
& matrix(1, 3) * point(3) + matrix(1, 4)
        temp(2) = matrix(2, 1) * point(1) + matrix(2, 2) * point(2) + 
& matrix(2, 3) * point(3) + matrix(2, 4)
        temp(3) = matrix(3, 1) * point(1) + matrix(3, 2) * point(2) + 
& matrix(3, 3) * point(3) + matrix(3, 4)
        temp(4) = matrix(4, 1) * point(1) + matrix(4, 2) * point(2) + 
& matrix(4, 3) * point(3) + matrix(4, 4)

        call ptkf_pt4topt3(temp, tpoint)
        
        RETURN
        END
        
        SUBROUTINE ptkf_transform(matrix, point, tpoint)
C /*
C ** \parambegin
C ** \param{REAL}{matrix(3, 3)}{3x3 matrix}{IN}
C ** \param{REAL}{point(2)}{2D point}{IN}
C ** \param{REAL}{tpoint(2)}{2D point}{OUT}
C ** \paramend
C ** \blurb{This function transforms the 2D point \pardesc{point} by the
C ** $3 \times 3$ matrix \pardesc{matrix} and returnes the result
C ** in {\tt tpoint}.}
C */
        REAL matrix(3,3), point(2), tpoint(2)

        tpoint(1) = matrix(1, 1) * point(1) + matrix(1, 2) * point(2) + 
& matrix(1, 3)
        tpoint(2) = matrix(2, 1) * point(1) + matrix(2, 2) * point(2) + 
& matrix(2, 3)
        
        RETURN
        END
        
        SUBROUTINE ptkf_matrixtomatrix3(mat, mat3)
C /*
C ** \parambegin
C ** \param{REAL}{mat(3, 3)}{3x3 matrix}{IN}
C ** \param{REAL}{mat3(4, 4)}{4x4 matrix}{OUT}
C ** \paramend
C ** \blurb{This function converts the $3 \times 3$ matrix \pardesc{mat}
C ** to the $4 \times 4$ matrix \pardesc{mat3}, as follows:
C ** $$  \left( \begin{array}{ccc}
C **          a & b & c\\
C **          d & e & f\\
C **          g & h & j
C **        \end{array} \right)
C ** 
C ** \rightarrow
C **
C ** \left( \begin{array}{cccc}
C **          a & b & 0 & c\\
C **          d & e & 0 & f\\
C **          0 & 0 & 1 & 0\\
C **          g & h & 0 & j
C **        \end{array} \right) 
C ** $$}
C */
        REAL mat(3,3), mat3(4,4)
        external ptk_matrixtomatrix3 !$PRAGMA C(ptk_matrixtomatrix3)

        call  ptk_matrixtomatrix3(mat, mat3)
        
        RETURN
        END
        
        SUBROUTINE ptkf_outputmatrix3(fileptr, matrix, string)
C /*
C ** \parambegin
C ** \param{INTEGER}{fileptr}{file pointer}{OUT}
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{IN}
C ** \param{CHARACTER*(*)}{string}{character string}{IN}
C ** \paramend
C ** \blurb{This function outputs the $4\times 4$ matrix
C **  \pardesc{matrix} and the message \pardesc{string}
C ** to the file specified by \pardesc{fileptr}.}
C */
        INTEGER fileptr
        REAL matrix(4,4)
        CHARACTER*(*) string
        external ptk_outputmatrix3 !$PRAGMA C(ptk_outputmatrix3)
        
        call  ptk_outputmatrix3(%val(fileptr), matrix, string)
        
        RETURN
        END
        
        SUBROUTINE ptkf_box3tobox3(box1, box2, preserve, operation,
& matrix, error)
C /*
C ** \parambegin
C ** \param{REAL}{box1(6)}{3D volume}{IN}
C ** \param{REAL}{box2(6)}{3D volume}{IN}
C ** \param{LOGICAL}{preserve}{preserve aspect ratio}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{OUT}
C ** \param{INTEGER}{error}{error code}{OUT}
C ** \paramend
C ** \blurb{This function computes a mapping from one 3D
C ** box to another -- a 3D window to 3D viewport
C ** transformation -- and concatenates this transformation
C ** with \pardesc{matrix} on the
C ** basis of \pardesc{operation}.If the parameters are invalid,
C **  \pardesc{error} is set to 
C ** -1. Otherwise, its value is 0.}
C */
        REAL box1(6), box2(6)
        LOGICAL preserve
        INTEGER operation
        REAL matrix(4,4)
        INTEGER error
        external ptk_box3tobox3 !$PRAGMA C(ptk_box3tobox3)

        call  ptk_box3tobox3(box1, box2, %val(preserve), 
& %val(operation), matrix, error)
        
        RETURN
        END
        
        SUBROUTINE ptkf_boxtobox(box1, box2, preserve, operation, 
& matrix, error)
C /*
C ** \parambegin
C ** \param{REAL}{box1(4)}{2D box}{IN}
C ** \param{REAL}{box2(4)}{2D box}{IN}
C ** \param{LOGICAL}{preserve}{preserve aspect ratio}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(3, 3)}{3x3 matrix}{OUT}
C ** \param{INTEGER}{error}{error code}{OUT}
C ** \paramend
C ** \blurb{This function computes
C **  a mapping from one 2D box to another -- a 2D window to 2D viewport
C ** transformation ---and concatenates this transformation
C ** it with \pardesc{matrix} on the
C ** basis of \pardesc{operation}.
C ** If the parameters are invalid, \pardesc{error} is set to 
C ** -1. Otherwise, its value is 0.}
C */
        REAL box1(4), box2(4)
        LOGICAL preserve
        INTEGER operation
        REAL matrix(3,3)
        INTEGER error
        external ptk_boxtobox  !$PRAGMA C(ptk_boxtobox)

        call ptk_boxtobox(box1, box2, %val(preserve), %val(operation), 
& matrix, error)
        
        RETURN
        END
        
        SUBROUTINE ptkf_accumulatetran3(fixed, shift, rotx, roty, rotz, 
& scale,operation, matrix)
C /*
C ** \parambegin
C ** \param{REAL}{fixed(3)}{origin}{IN}
C ** \param{REAL}{shift(3)}{shift factor}{IN}
C ** \param{REAL}{rotx}{x rotation}{IN}
C ** \param{REAL}{rotx}{y rotation}{IN}
C ** \param{REAL}{rotx}{z rotation}{IN}
C ** \param{REAL}{scale(3)}{scale factor}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{OUT}
C ** \paramend
C ** \blurb{This function computes the specified 3D
C ** shift, scale and rotate transformations, in the order 
C ** scale, rotate, shift, 
C ** and then concatenates the resulting transformation
C ** with the specified matrix on the basis of \pardesc{operation}.}
C */
        REAL fixed(3), shift(3), rotx, roty, rotz, scale(3)
        INTEGER operation
        REAL matrix(4,4)
        REAL*8 dprotx, dproty, dprotz
        external ptk_accumulatetran3 !$PRAGMA C(ptk_accumulatetran3)
        
        dprotx = rotx
        dproty = roty
        dprotz = rotz
        call  ptk_accumulatetran3(fixed, shift, %val(dprotx), 
& %val(dproty), %val(dprotz), scale, %val(operation), matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_accumulatetran(fixed, shift, rot, scale, 
& operation, matrix)
C /*
C ** \parambegin
C ** \param{REAL}{fixed(2)}{origin}{IN}
C ** \param{REAL}{shift(2)}{shift factor}{IN}
C ** \param{REAL}{rotx}{x rotation}{IN}
C ** \param{REAL}{scale(2)}{scale factor}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(3, 3)}{3x3 matrix}{OUT}
C ** \paramend
C ** \blurb{This function computes the specified 2D
C ** shift, scale and rotate transformations, in the order 
C ** scale, rotate, shift, 
C ** and then concatenates the resulting transformation
C ** with the specified matrix on the basis of \pardesc{operation}.}
C */
        REAL fixed(2), shift(2), rot, scale(2)
        INTEGER operation
        REAL matrix(3,3)
        REAL*8 dprot
        external ptk_accumulatetran !$PRAGMA C(ptk_accumulatetran)
        
        dprot = rot
        call  ptk_accumulatetran(fixed, shift, %val(dprot), scale, 
& %val(operation), matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_evalvieworientation3(viewrefpoint, 
& viewplanenormal, viewupvector, operation, matrix, error)
C /*
C ** \parambegin
C ** \param{REAL}{viewrefpoint(3)}{view reference point}{IN}
C ** \param{REAL}{viewplanenormal(3)}{view plane normal}{IN}
C ** \param{REAL}{viewupvector(3)}{view up vector}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{OUT}
C ** \param{INTEGER}{error}{error code}{OUT}
C ** \paramend
C ** \blurb{This function computes 
C ** a 3D PHIGS view orientation matrix on the basis of
C ** a specified view reference point (\pardesc{viewrefpoint}), a
C ** view plane normal (\pardesc{viewplanenormal}) and a view up vector
C ** (\pardesc{viewupvector}). If the function succeeds,
C **  \pardesc{error} is set to 
C **  0. Otherwise,
C ** \pardesc{error} is 61 if the view plane normal is null,
C ** 63 if the view up vector is null,
C ** and 58 if the cross product of the view up vector
C ** and the view plane normal is null.}
C */
        REAL viewrefpoint(3), viewplanenormal(3), viewupvector(3)
        INTEGER operation
        REAL matrix(4,4)
        INTEGER error
        external ptk_evalvieworientation3 
& !$PRAGMA C(ptk_evalvieworientation3)
        
        call  ptk_evalvieworientation3(viewrefpoint, viewplanenormal, 
& viewupvector, %val(operation), matrix, error)
        
        RETURN
        END
        
        SUBROUTINE ptkf_evalvieworientation(viewrefpoint, viewupvector, 
& operation, matrix, error)
C /*
C ** \parambegin
C ** \param{REAL}{viewrefpoint(2)}{view reference point}{IN}
C ** \param{REAL}{viewupvector(2)}{view up vector}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(3, 3)}{3x3 matrix}{OUT}
C ** \param{INTEGER}{error}{error code}{OUT}
C ** \paramend
C ** \blurb{This function computes 
C ** a 2D PHIGS view orientation matrix on the basis of
C ** a specified view reference point (\pardesc{viewrefpoint})
C ** and a view up vector
C ** (\pardesc{viewupvector}).  If the function succeeds,
C ** \pardesc{error} is set to 0. Otherwise, 
C ** \pardesc{error} is 63 if the view up vector is null.}
C */
        REAL viewrefpoint(2), viewupvector(2)
        INTEGER operation
        REAL matrix(3,3)
        INTEGER error
        external ptk_evalvieworientation 
& !$PRAGMA C(ptk_evalvieworientation)

        call  ptk_evalvieworientation(viewrefpoint, viewupvector, 
& %val(operation), matrix, error)
        
        RETURN
        END
        
        SUBROUTINE ptkf_evalviewmapping3(wlimits, vlimits, viewtype, 
& refpoint,vplanedist, operation, matrix, error)
C /*
C ** \parambegin
C ** \param{REAL}{wlimits(6)}{window limits}{IN}
C ** \param{REAL}{vlimits(6)}{viewport limits}{IN}
C ** \param{INTEGER}{viewtype}{projection type}{IN}
C ** \param{REAL}{refpoint(3)}{projection reference point}{IN}
C ** \param{REAL}{vplanedist}{view plane distance}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{OUT}
C ** \param{INTEGER}{error}{error code}{OUT}
C ** \paramend
C ** \blurb{This function evaluates a 3D PHIGS view mapping matrix.
C ** If the function succeeds,
C ** \pardesc{error} is set to 0. Otherwise, 
C ** \pardesc{error} is 
C ** 329 if the window limits are not valid,
C ** 336 if the back plane is in front of front plane,
C ** 330 if the viewport limits are not valid,
C ** 335 if the projection reference point is on the view plane,
C ** and 340 if the projection reference point is between front and back planes.}
C */
        REAL wlimits(6), vlimits(6)
        INTEGER viewtype
        REAL refpoint(3), vplanedist
        INTEGER operation
        REAL matrix(4,4)
        INTEGER error
        REAL*8 dpvplanedist
        external ptk_evalviewmapping3 !$PRAGMA C(ptk_evalviewmapping3)

        dpvplanedist = vplanedist
        call  ptk_evalviewmapping3(wlimits, vlimits, %val(viewtype), 
& refpoint,%val(dpvplanedist), %val(operation), matrix, error)
        
        RETURN
        END
        
        SUBROUTINE ptkf_evalviewmapping(wlimits, vlimits, operation, 
& matrix, error)
C /*
C ** \parambegin
C ** \param{REAL}{wlimits(4)}{window limits}{IN}
C ** \param{REAL}{vlimits(4)}{viewport limits}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(3, 3)}{3x3 matrix}{OUT}
C ** \param{INTEGER}{error}{error code}{OUT}
C ** \paramend
C ** \blurb{This function evaluates a 2d PHIGS view mapping matrix.
C ** If the function succeeds,
C ** \pardesc{error} is set to \pardesc{ptkcpcok}. Otherwise, 
C ** \pardesc{error} is
C ** 329 if the window limits are not valid,
C ** and 330 if the viewport limits are not valid.}
C */
        REAL wlimits(4), vlimits(4)
        INTEGER operation
        REAL matrix(3,3)
        INTEGER error
        external ptk_evalviewmapping !$PRAGMA C(ptk_evalviewmapping)

        call  ptk_evalviewmapping(wlimits, vlimits, %val(operation), 
& matrix, error)
        
        RETURN
        END
        
        SUBROUTINE ptkf_stackmatrix3(matrix)
C /*
C ** \parambegin
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{IN}
C ** \paramend
C ** \blurb{This function pushes
C ** the $4 \times 4$ matrix \pardesc{matrix} onto the transformation stack.}
C */
        REAL matrix(4,4)
        external ptk_stackmatrix3 !$PRAGMA C(ptk_stackmatrix3)

        call  ptk_stackmatrix3(matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_stackmatrix(matrix)
C /*
C ** \parambegin
C ** \param{REAL}{matrix(3, 3)}{3x3 matrix}{IN}
C ** \paramend
C ** \blurb{This function pushes the $3 \times 3$ matrix \pardesc{matrix}
C **  onto the transformation stack.}
C */
        REAL matrix(3,3)
        external ptk_stackmatrix !$PRAGMA C(ptk_stackmatrix)

        call  ptk_stackmatrix(matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_unstackmatrix3(matrix)
C /*
C ** \parambegin
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{OUT}
C ** \paramend
C ** \blurb{This function
C ** pops a $4 \times 4$  matrix 
C ** from the transformation stack and returns it in
C ** \pardesc{matrix}.}
C */
        REAL matrix(4,4)
        external ptk_unstackmatrix3 !$PRAGMA C(ptk_unstackmatrix3)

        call  ptk_unstackmatrix3(matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_unstackmatrix(matrix)
C /*
C ** \parambegin
C ** \param{REAL}{matrix(3, 3)}{3x3 matrix}{OUT}
C ** \paramend
C ** \blurb{This function pops a $3 \times 3$ matrix
C ** from the transformation stack and returns it in
C ** \pardesc{matrix}.}
C */
        REAL matrix(3,3)
        external ptk_unstackmatrix !$PRAGMA C(ptk_unstackmatrix)

        call  ptk_unstackmatrix(matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_examinestackmatrix3(matrix)
C /*
C ** \parambegin
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{IN}
C ** \paramend
C ** \blurb{This function returns the top entry on the transformation stack.
C ** The stack is not disturbed.}
C */
        REAL matrix(4,4)
        external ptk_examinestackmatrix3 
& !$PRAGMA C(ptk_examinestackmatrix3)

        call  ptk_examinestackmatrix3(matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_examinestackmatrix(matrix)
C /*
C ** \parambegin
C ** \param{REAL}{matrix(3, 3)}{3x3 matrix}{IN}
C ** \paramend
C ** \blurb{This function  returns the top entry on the transformation stack.
C ** The stack is not disturbed.}
C */
        REAL matrix(3,3)
        external ptk_examinestackmatrix 
& !$PRAGMA C(ptk_examinestackmatrix)

        call  ptk_examinestackmatrix(matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_3ptto3pt(p1, p2, p3, q1, q2, q3, operation, 
& matrix, error)
C /*
C ** \parambegin
C ** \param{REAL}{p1(3)}{3D point}{IN}
C ** \param{REAL}{p2(3)}{3D point}{IN}
C ** \param{REAL}{p3(3)}{3D point}{IN}
C ** \param{REAL}{q1(3)}{3D point}{IN}
C ** \param{REAL}{q2(3)}{3D point}{IN}
C ** \param{REAL}{q3(3)}{3D point}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{OUT}
C ** \param{INTEGER}{error}{error code}{OUT}
C ** \paramend
C ** \blurb{This function returns the 3 point to 3 point transformation as
C ** described in \cite{mort:geom}, pages 353--355.
C ** The transformation has the following properties:
C ** \pardesc{p1} is transformed onto \pardesc{q1};
C ** the vector \pardesc{(p2-p1)} is transformed to be parallel to the vector 
C ** \pardesc{(q2-q1)};
C ** the plane containing the three points \pardesc{p1, p2, p3} is
C ** transformed into  the plane containing \pardesc{q1, q2, q3}.
C ** The transformation is concatenated with the $4 \times 4$ matrix
C ** \pardesc{matrix} on the basis of \pardesc{operation}.
C ** If the parameters are invalid, \pardesc{error} is set to 
C ** -1. Otherwise, its value is \pardesc{ptkcpcok}.}
C */
        REAL p1(3), p2(3), p3(3), q1(3), q2(3), q3(3)
        INTEGER operation
        REAL matrix(4,4)
        INTEGER error
        external ptk_3ptto3pt  !$PRAGMA C(ptk_3ptto3pt)

        call  ptk_3ptto3pt(p1, p2, p3, q1, q2, q3, %val(operation), 
& matrix, error)
        
        RETURN
        END
        
        SUBROUTINE ptkf_0to3pt(origin, xdirn, ydirn, operation, matrix)
C /*
C ** \parambegin
C ** \param{REAL}{origin(3)}{origin of axes}{IN}
C ** \param{REAL}{xdirn(3)}{x direction}{IN}
C ** \param{REAL}{y dirn(3)}{y direction}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{OUT}
C ** \paramend
C ** \blurb{This function computes an object transformation which
C ** maps unit vectors 
C ** along the $x$, $y$ and $z$ axes onto unit vectors along the
C ** corresponding axes 
C ** of the new coordinate system.}
C */
        REAL origin(3), xdirn(3), ydirn(3)
        INTEGER operation
        REAL matrix(4,4)
        external ptk_0to3pt !$PRAGMA C(ptk_0to3pt)
        
        call  ptk_0to3pt(origin, xdirn, ydirn, %val(operation), matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_oto3pt(origin, xdirn, ydirn, operation, matrix)
C /*
C ** \parambegin
C ** \param{REAL}{origin(3)}{origin of axes}{IN}
C ** \param{REAL}{xdirn(3)}{x direction}{IN}
C ** \param{REAL}{y dirn(3)}{y direction}{IN}
C ** \param{INTEGER}{operation}{concatenation operation}{IN}
C ** \param{REAL}{matrix(4, 4)}{4x4 matrix}{OUT}
C ** \paramend
C ** \blurb{This function performs the same operation as
C **  \pardesc{ptk\_0to3pt}, except the name has an \pardesc{o}\ 
C ** (oh) instead of \pardesc{0}\ (zero). This function is provided for members
C ** of the Fumbly Fingers Club.}
C */
        REAL origin(3), xdirn(3), ydirn(3)
        INTEGER operation
        REAL matrix(4,4)
        external ptk_oto3pt !$PRAGMA C(ptk_oto3pt)        

        call  ptk_oto3pt(origin, xdirn, ydirn, %val(operation), matrix)
        
        RETURN
        END
        
        SUBROUTINE ptkf_invertmatrix3(a, ainverse, error)
C /*
C ** \parambegin
C ** \param{REAL}{a(4, 4)}{4x4 matrix}{IN}
C ** \param{REAL}{ainverse(4, 4)}{4x4 matrix}{OUT}
C ** \param{INTEGER}{error}{error code}{OUT}
C ** \paramend
C ** \blurb{This function computes the inverse of the $4 \times 4$
C ** matrix \pardesc{a}, 
C ** returning the result in \pardesc{ainverse}. 
C ** If matrix \pardesc{a} is singular, then 
C ** \pardesc{error} is set to $-1$, and \pardesc{ainverse} is undefined,
C ** otherwise \pardesc{error} is set to 0.}
C */
        REAL a(4,4), ainverse(4,4)
        INTEGER error
        external ptk_invertmatrix3 !$PRAGMA C(ptk_invertmatrix3)
        
        call  ptk_invertmatrix3(a, ainverse, error)
        
        RETURN
        END
        
        SUBROUTINE ptkf_invertmatrix(a, ainverse, error)
C /*
C ** \parambegin
C ** \param{REAL}{a(3, 3)}{3x3 matrix}{IN}
C ** \param{REAL}{ainverse(3, 3)}{3x3 matrix}{OUT}
C ** \param{INTEGER}{error}{error code}{OUT}
C ** \paramend
C ** \blurb{This function computes the inverse of the $3 \times 3$
C ** matrix \pardesc{a}, 
C ** returning the result in \pardesc{ainverse}. 
C ** If matrix \pardesc{a} is singular, then 
C ** \pardesc{error} is set to $-1$, and \pardesc{ainverse} is undefined,
C ** otherwise \pardesc{error} is set to 0.}
C */
        REAL a(4,4), ainverse(4,4)
        INTEGER error
        external ptk_invertmatrix !$PRAGMA C(ptk_invertmatrix)
        
        call  ptk_invertmatrix(a, ainverse, error)
        
        RETURN
        END
        
C       end of tran.f 
